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FACTORIZATION WITH GENUS 2 CURVES 

ROMAIN COSSET 

Abstract. The elliptic curve method (ECM) is one of the best factorization 
methods available. It is possible to use hyperelliptic curves instead of elliptic 
curves but it is in theory slower. We use special hyperelliptic curves and Rum- 
mer surfaces to reduce the complexity of the algorithm. Our implementation 
GMP-HECM is faster than GMP-ECM for factoring big numbers. 



Introduction 

The elliptic curve method (ECM) introduced in 1985 by H. W. Lenstra, Jr. [6] 
plays an important role in factoring integers. ECM is used to find "medium sized" 
(up to 60 digits) prime factors of "random" numbers. It is also used by other fac- 
toring algorithms like sieving methods for cofactoring [5]. ECM is a generalization 
of Pollard's p — 1 algorithm: instead of working in F*, we work in the group of 
points on an elliptic curve. We generalize it by using hyperelliptic curves of genus 2 
instead of elliptic curves of genus 1. 

At first sight the "hyperelliptic curve method" (HECM) seems slower that ECM 
because of two reasons: first the arithmetic of hyperelliptic curves is slower com- 
pared to the arithmetic of elliptic curves, and secondly its probability of success is 
smaller than that of ECM. 

Indeed this probability depends on the probability of the cardinality of the Jacobian 
of the curve modulo a prime factor p of n being smooth. But the probability of a 
number being smooth decreases with its size. By Weil's theorem, the cardinality 
of a Jacobian of a genus 2 curve on F p is around p 2 whereas the cardinality of an 
elliptic curve over the same field is around p. This seems like a major deterrent 
to using HECM. To avoid it we use special hyperelliptic curves whose Jacobians 
are isogenous to the product of two elliptic curves: they are called decomposable. 
One run of HECM with this kind of curves is equivalent to two simultaneous runs 
of ECM, thus the probability of success of HECM with decomposable hyperelliptic 
curves is comparable to that of two ECM. 

The complexity of stage 1 is dominated by the cost of the arithmetic on the curves. 
In genus 1, the quickest arithmetic is obtained by using Montgomery formulae [8] or 
the new Edwards curves [1]. In genus 2, the best known formula? use Kummer sur- 
faces [3]. The Kummer surface is a variety obtained by identifying opposite points 
of the Jacobian. However, not all genus 2 hyperelliptic curves map to a rational 
Kummer surface. 
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The curves used in HECM arc built as reduction modulo n of curves denned over 
Q which are selected for satisfying different constraints. We need a large number 
of possible trial curves in order to factor n, hence the requirement that there exists 
an infinite family of curves over Q. This requires some work; finally we obtain a 
parametrization of a subfamily of hyperelliptic curves: the parameters live on an 
elliptic curve over a function field. 

We have implemented our algorithm by using many functions of GMP-ECM a 
free ECM program. During the computation, there are many multiplications by 
parameters of the hyperelliptic curves so the use of small parameters (i.e. which 
fit into one machine word) makes these multiplications negligible before the cost 
of full-length modular multiplications. This makes our software faster than GMP- 
ECM. 

In Section 1 we give some facts about ECM, decomposable genus 2 hyperelliptic 
curves and Kummer surfaces. Details of the parametrization are given in Section 
2 while Section 3 focuses on the HECM algorithm. In Section 4 we describe our 
implementation and give some numerical results. 

In the following, n denotes the number to be factored. 

1. BACKGROUND 

1.1. ECM. Following a common notation abuse in factoring algorithms, we work 
over Z/nZ as if it were a field. The only operation that might fail is "field" inversion 
which is calculated using the Euclidean algorithm. If an inversion fails, we find a 
factor of n. 

The ECM method starts by choosing a random elliptic curve £ over Z/nZ and 
a point P on it. Let He a positive integer, we compute Q = [k]P and we hope Q 
to be the zero of the curve modulo a prime p dividing n but not modulo n. If this 
is the case, then one division fails during the computation and we find the factor p. 
This happens if the cardinality \£ (F p ) | divides k. There are two different phases 
(also called stages or steps) in ECM: in phase 1 we hope that the cardinality is 
i?i-smooth (i.e. multiple of prime powers less than Si) for some bound B\. In 
phase 2 we try to cover the "close-miss" case where the cardinality is a product of 
a i?i-smooth part by a prime cofactor. If this fails, we take another elliptic curve. 
ECM is a probabilistic algorithm whose complexity is 0(L(p)^ 2+ °^ M(log(n))) 
with L(p) = e V lo s(p) i°g(i°g(p)) an d M(log(n)) is the complexity of multiplications 
modulo n. The complexity of ECM is dominated by the size of the smallest factor 
p of n rather than the size of the number n to be factored. Note however that ECM 
does not always find the smallest factor. 

In stage 1 we hope that the cardinality of the elliptic curve is Si-smooth. We 
take k = \\^ Bl ^ Xo ^ B ^/ Xo ^ = 1cm (1,2,..., B x ) so that all Bi-smooth numbers 
divide k. The main cost of stage 1 is the arithmetic; different methods are used to 
reduce this cost. We focus on ECM with Montgomery formulae [8] since it shares 
several traits with the HECM method which is the subject of this work. Let £ 
be an elliptic curve in Montgomery's form: y 2 = x 3 + ax 2 + x. We use projective 
coordinates (x : y : z) to avoid divisions which are expensive. Moreover we disregard 
the y-coordinate which means that we identify a point with its opposite. It is still 
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possible to "double" a point (i.e. computing ±2P) and if the difference ± (P — Q) 
between two points is known then we can compute their sum ±(P + Q): this 
operation is called pseudo-addition. To compute the multiple of a point we need 
to find chains of doubling and pseudo-additions which is a special case of addition 
chains called "Lucas Chains" [7]. For instance 

1^2^3^4^7^10^17 

is a chain for 17. One way to find such chains is to note that if we know [n]P and 
[n + 1] P then we can compute [2n] P, [2n + 1] P and [2n + 2] P. Hence we have 
binary chains: at each step we choose the point to double according to the binary 
expansion of the multiplier. For instance the following chain is the binary chain 
for 17: 

1^2^3^4^5^8^9^17. 

Binary chains are not the shortest chain of doubling and pseudo-addition. Mont- 
gomery's PRAC algorithm finds short Lucas chains [7]. 

A description of stage 2 falls out of scope of the present article since for reasons 
which will be explained later, stage 1 is the only phase done by HECM. 

1.2. Decomposable curves. 

Definition 1.1. A genus 2 hyperelliptic curve C is said to be decomposable if there 
exists an isogeny <p from Jac(C) to the product of two elliptic curves Ec 

<p : Jac (C) — > £i x £ 2 . 

If the kernel of the isogeny is Z/kZ x Z/kZ we say that the curve is (k,k)- 
decomposable. 

There are characterizations for a curve to be (fc, ^-decomposable [11]. We will 
only look at (2, 2)-decomposable curves since the conditions are the simplest. 

Theorem 1.2. Let C be a hyperelliptic curve of genus 2 given by the equation 

y 2 = x (x — 1) (x — A) (x — /i) (x — v) 

where \, [i, v are the Rosenhain invariants of the curve. C is (2, 2) -decomposable 
if and only if its Rosenhain invariants are linked by 



The proof is given in [4] . The two underlying elliptic curves have equation 
XV 2 = (x - 1) (x - x%) (x - xf) 

with 

q = ±vV(m - v )i x 2 = , x 3 = — — — , x = -quip - 1) • 

a* q l a* ~r q 

Over non algebraically closed fields K, the Rosenhain invariants may be non ra- 
tional. Thus a (2, 2)-decomposable curve can have a different equation. If C is a 
(2, 2)-decomposable curve in Rosenhain form then the two underlying elliptic curves 
are rational if and only if \i {\i — v) is a square in K. 
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The morphisms between the hyperelliptic curves and the elliptic curves are defined 
on K (q) . They are given by 

1.3. Kummer surfaces. Gaudry proposes in [3] to use theta functions to perform 
the arithmetic in the Jacobian of some genus 2 curves. Working with theta functions 
allows the use of the numerous classical formulae that were found using complex 
analysis: the reader is referred to Mumford's books [9, 10] for details. In fact all 
the formulae we need are algebraic and can be used on finite fields (of characteristic 
different from 2). 

Let n be a matrix in the Siegel half space of dimension 2. The theta functions are 
functions from C 2 / (Z 2 + f2Z 2 ) to C. There are 16 theta functions with half-integer 
characteristic: 10 are even and 6 are odd. The scalars obtained by evaluating them 
in z = (0, 0) are called theta constants. The main difference with Gaudry's article 
is that we use the squares of the coordinates. We are principally interested in 
only four squares of theta functions and the corresponding square theta constants 
(a : (i : 7 : 5). 

A Kummer surface JC^-p.-y-s) can be defined by the choice of four "theta constants" 
(a : (3 : 7 : 5) in P 3 (C). We write (X : Y : Z : T) for the coordinates of a point on 
fc(a:f3:-y.8)- The equation of the Kummer surface is: 

4E' 2 af3~fSXYZT = ((X 2 + Y 2 + Z 2 + T 2 ) - 

-F (XT + YZ) - G (XZ + YT) - H (XY + ZT)f 

where E' , F, G, H are constants that can be calculated from (a : (3 : 7 : 5) by the 
formulae given in the appendix. 

Let 4 be the ratio of two theta constants which can be calculated from the theta 



constants (a : (3 : 7 : 5). Let 

76 ae 

X: =—y ^'"w v '- = w 

C : y 2 = x (x — 1) (x — A) (x — ^1) (x — v) . 
Then Jac(C) / {±1} is isomorphic to JC^p.^,^ over C: 

■0: Jac(C)/{±l} — > K {a:l3 ,^. s) 

Over non algebraically closed fields, let C be the quadratic twist of C. Then 
Jac(C) /{±1} and Jac(C)/ {±1} are included in JC^ a .fj. T sy. 

JC {a :p: T .s) - V (Jac (C) I {±1}) |J i> (jac (c) / {±1}) 

with the 2-torsion points shared by the two Jacobians. Usually, we use the Mumford 
coordinates (u, v) on Jac (C), from a point P in K( a: p :T .g), it is possible to calculate 
u and v 2 of its image by the morphism [3, 12]. Of course v is defined up to its sign 
since we can't distinguish between one point and its opposite. This generalizes the 
idea of Montgomery's formulae for elliptic curves. 

The arithmetic on Jac(C) transports to arithmetic on /C( Q:( g :7:( 5). Given a point 
P = ip (D) on the Kummer surface it is possible to double it (i.e. to compute 
tp ([2]D)): see algorithm 1. Given two points on a Kummer surface P = tp (D) and 
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Algorithm 1 Doubling on a Kummer surface 



Input: a point P = (X : Y : Z : T) on K 


(a:f3: T .6)- 








Output: the point [2]P = (X 2 :Y 2 : Z 2 : 


T2). 








X' = (X + Y + Z + T) 2 i, 


Y' 


= (X + Y - 


■Z-T) 2 ± 


! 


Z' = {X-Y + Z-Tf±, 


r 


= [X-Y- 


-Z + T) 2 ± 




x 2 = (x' + Y' + z' + rf^, 


Y 2 


= (X' + Y' 


_ Z I _ T /)2 


1 

/3' 


Z 2 = (X' -Y' + Z' -T') 2 i, 


T 2 


= (X' - Y' 


-Z' + T'f 


1 

.5 



Algorithm 2 Pseudo-addition on a Kummer surface 

Input: two points P = (X : Y : Z : T) and Q = (X : Y : Z : T) on /C (ct:/3:7:5) , and 

the point R= (X : Y : Z : f) equal to P - Q such that XYZT ^ 0. 
Output: the point P + Q = (x : y : z : t). 



X' 


= (X + Y + Z + 


T)(X 


+ Y + 


Z + T)± 
Z~T)% 
Z~T)l, 
Z + T)l 


Y' 


= (X + Y-Z- 


T)(X 


+ Y- 


Z' 


= (X-Y + Z- 


T)(X 


-Y + 


V 


= (X-Y-Z + 


T)(X 

+ rf 

-T'f 
-T'f 


-Y- 


X 


= (X' + Y' + Z' 


i 

X ' 




y 


= (X' + Y'-Z' 


i 

Y> 




z 


= (X' -Y' + Z' 


1 

Z> 




t 


= (X' -Y'-Z' 


+ T'f 


1 

T 





Q = ip(D'), we can't add them since we don't know whether we must compute 
tp (D + D') or ip (D — D') . However if one of these quantities is known then it is 
possible to compute the other, this is called a pseudo-addition (algorithm 2). We 
note [2]P and P + Q these two operations. Note that these formulae do not hold 
when one coordinate is zero. 

We work with projective coordinates so divisions can be replaced by multipli- 
cations. Moreover the constants 1/a, 1/(3, I/7, l/S, 1/A, 1/B, l/C and 1/D can 
be precomputed. The cost of pseudo-addition is 4 divisions, 4 multiplications, 4 
squares and 4 multiplications by constants (thereafter denoted 41 + AM + 45 + 4rf) 
or if the divisions are replaced by multiplications the cost becomes 1AM + 45 + Ad. 
The cost of doubling is 85 + 8d. Note that by using the properties of projective 
coordinates, it is possible to save some multiplications [3]. 

For the multiplication algorithm we want to avoid divisions since they are costly. 
This imposes that the inverses of the coordinates of P — Q are always known when 
we want to pseudo-add P and Q. Thus it is impossible to use "short" Lucas chains 
and the PRAC algorithm, we must use binary chains where P — Q is always the 
initial point (algorithm 3). For each bit of the multiplier k, we do one doubling 
and one pseudo-addition on the Kummer surface. The total cost of a multiplication 
by k is 7M + 95 + 9d per bit of k since some computations can be shared between 
the two operations (it is 12M + 45 + 16g? if we don't share them). 

Remark 1.3. The coordinates of the initial point P affect the cost of computation. 
For instance if two coordinates of P are equal or opposite then we gain, for each 
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Algorithm 3 Multiplication algorithm 

Input: a point P on )C( a: p :i: s) and an integer k. 
Output: the point [k]P. 

if k = 2 then 
return [2}P 

else 

Let k = J2\=o be the binary expansion of k with fc = 1 the most 

significant bit. 

P m =P;P P = [2]P; 

for i from 2 to / do 

Q = P p + P m {note that we have P p — P rn = P} 

if fcj = f then 

Pp — [2]-Pp! = Q\ 

else {fci = 0} 

Pm [2]-Fm? Pp Q: 

end if 
end for 
return P m 
end if 



bit of the multiplier, one multiplication and another one which in fact is a square. 
This is a speed up of 5%. 

2. Parametrization 

In this Section, we find suitable curves for HECM. We obtain curves over Q 
having the desired properties and consider their reduction modulo n. We want to 
have an infinite number of curves over Q so that we have enough curves over Z/nZ. 
Note that throughout the construction, we should never compute a square root be- 
cause of two reasons: first this square root could be in an extension of Q, thus when 
reducing modulo n wc could be working in an extension of Z/nZ and the arithmetic 
would be slower. Moreover, computing square roots modulo the composite integer 
n is equivalent to find factors of n. To sum up, we require that all constants are 
defined by rational functions in several variables. 

2.1. Parametrization of the hyperelliptic curves. 

Lemma 2.1. Let C be a genus 2 curve over Q with equation 

C : XV 2 = x (x — 1) (x — X) (x — /«) (x — v) . 

C can be used in HECM (i.e. C is (2, 2) -decomposable with rational underlying 
elliptic curves and the rational Kummer surface) if and only if 

A = fi- — — , //(// — !/)=□, Xfif = D 
1-fj, 

where □ means that the quantity must be a square. 

Proof. The first two conditions mean that C is (2, 2)-decomposable with rational 
underlying elliptic curves. The Kummer surface is rational if and only if the squares 
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(a : [3 : 7 : 6) of the four theta constants are rational. They are linked with the 
Roscnhain invariants by: 

0:7 7e ae 

where e and </> are two other squares of theta constants whose ratio is rational when 
(a : (3 : 7 : S) are rational. Note that this implies that the Rosenhain invariants 
must also be rational so it gives a justification to the choice of the equation of the 
hyperelliptic curve. From these three equations, we obtain: 



a s/Xfjbv 7 



(3 yU 8 v 

So the product \[iv must be a square in Q. 

When we write the equation A = ^jE^ m terms of theta constants (use the formulae 
from [4] and [12]), it yields to a 2 = S 2 . Moreover (3, 7, S are functions of A, fi, v 
and a and thus are rational. This proves that the conditions arc sufficient. □ 

Remark 2.2. The Rosenhain invariants (A,^, v) are defined up to the action of the 
group PGL(2, 5) see [4]. Our choice here is one that leads to an equality between 
two of the first four theta constants. With a different choice of order, we would 
have had another square root to handle. Moreover the fact that a 2 = 5 2 is a main 
advantage for the arithmetic: it saves two scalar multiplications per multiplier's bit 
when computing the multiple of a point. 

Since (a : (3 : 7 : 5) live in P 3 (Q), we can take a = 1. We choose also 5 = a = 1. 
The choice 5 = — a = — 1 leads to an isomorphic Kummer surface. 
Take fi = 1 — '^^'^ with seQso that \[iv — [i 2 s 2 . The second equation becomes 

M (M - v) = {-v + v 2 + s 2 ) {u -l){v- s 2 ) = □. 

Assume that is a square u 2 (then v = s \zS ) so ^ nat tne equation above 

rewrites as: 

1 + (-3/s 2 + l/s> 2 + u A /s 2 = □. 

Let v 2 be the square, the point (u,v) lies on an elliptic curve over Q(s). This 
elliptic curve is in the Jacobi model [2], is of rank 1 and we have a non torsion 
point P = (l, 1 — -4-) on it. 

Theorem 2.3. A subfamily o/(2, 2) -decomposable hyperelliptic genus 2 curves with 
rational underlying elliptic curves and with rational Kummer surfaces is given by 
the following parametrization: 

C : XV 2 = x{x - 1) (x - A) (x - n) (x - v) , 

1 — v — s 2 — u 2 

X = l l -i 1 M=l o ' v= ~ l 2 

1 — /i s z 1 — u 

where (u, v) lies on the following elliptic curve 

1 + (-3/s 2 + l/s> 2 + u^/s 2 = v 2 . 

A non-torsion point on this curve is (1,1 — -^). The parameters s, \, u, v are 
rational numbers that must verify the following conditions. 
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Genericity condition 2.4. With the notations of the theorem, the curve C is of 
genus 2 if and only if 0, 1, oo, A, \i and v are distinct and \ 1S not zero. This is 
equivalent to 

X^O, s^0,±l, u96O.il, v^O, 
s ± ±u, s 2 - 2u 2 + u V 0. 

In particular, condition 2.4 implies that we can't use the point (1, 1 — p-) on 
the Jacobi curve. We must take a multiple of this point: for instance its double 
(2,1 + Jr). 

Genericity condition 2.5. For the arithmetic on the Kummer surface, all the 
theta constants must be non zero. The parameters must satisfy condition 2.4 and 

s £ ±u 2 . 

The parameter \ determines whether we are on the curve or on its quadratic 
twist. However it is not chosen during the paramctrization of the curve: since the 
curve and its twist are both included in the same Kummer surface, the choice of a 
point on it determines whether we work on the curve or on its twist. In practice, 
the y-coordinate of the point is never used, so x is never computed. 

We could have taken = lu 2 (here 1 = 1) but that would have led to more 
complicated equations. For instance with I = —1, we found an elliptic curve with 
no point on it. In this case (7 = —1), we had to assume that s was on a conic for 
the curve to become of rank 1. 

Remark 2.6. ECM supposes that the cardinality of elliptic curves on ¥ p behaves 
like a random integer of size around p with some additional divisibility conditions. 
In HECM, we work simultaneously with two elliptic curves. Their cardinalities 
are not independent: points of 4-torsion cannot exist on only one curve. How- 
ever, experimentally, it seems that this is the only noticeable link between their 
cardinalities. 

2.2. Finding a point on the Kummer surface. In ECM we need a generic 
point (i.e. a non torsion point) in the group. Contrary to the Brcnt-Suyama 
parametrization for elliptic curves, it is possible here to find a point after having 
chosen a hyperelliptic curve. However, the initial point P must be a non torsion 
point of the two underlying elliptic curves over Q and not only of the Jacobian of 
the hyperelliptic curve. Remember that we don't want to take square roots or to 
work in an extension field so we can't take three random coordinates and solve the 
equation of the Kummer surface to find the last one. 

The first method to find a point is to take one coordinate to be zero. Then the 
three other lie on a conic. Since conies are birationally equivalent to P 1 , we have 
an infinite number of points on the Kummer surface. In general, these points are 
not torsion points. Note however that these points cannot be used directly for the 
multiplication algorithm: we must double them until we find a point with no zero 
coordinate; in general, one doubling is enough. 

There is a better solution to find points: we look for points with two equal or 
opposite coordinates. Such points have the advantage of saving multiplications 
(see remark 1.3): the cost of the multiplication algorithm becomes 12M + 10S 
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per bit of the multiplier. It turns out that not all choices of coordinate pairs are 
possible: some choices lead only to torsion points or to impossible equations like 
— 1 = □. The choice Y = — X yields to the equation of an elliptic curve with 
many points on it. These points lead to points on the Kummer surface which are 
of infinite order (over Q). For instance we use the following point: 

2/2 n \ / 4 o 2, 2,3 2,\ 

u [u — l)[su —3su +u +s — s +sj y 
(s — u 2 )s 4 v 2 1 

T= (''-"'A" 3 - 1 ), Z = l. 

S V 

3. Algorithm and ameliorations 

3.1. HECM. In ECM, the big product k = 7r[1 ° s(Bl)/1 ° s(7r)1 is not im- 

puted as such [13]. Instead, for each prime ir we compute I — [log (Si) /log(7r)] 
and then we do / times Q — [tt]Q and go to the next prime. This avoids the 
computation of a big integer product and PRAC finds shorter chains if it works 
with prime multipliers [7] . If we work with Kummer surfaces we can't use PRAC 
because changing the initial point is costly (we would have divisions). Moreover 
the advantage of choosing a good initial point is lost if we don't use binary Lucas 
chains. Therefore in HECM we begin by computing k = U^<b 1 ttI 10 ^ 51 )/ 106 ^. 
In theory this computation is costly but it practice it is negligible: we use product 
tree and fast multiplication. Moreover this computation is shared for many curves. 

We hope to encounter the zero of one of the underlying elliptic curves so it is 
important to have explicit morphisms between the Kummer surface and the two 
underlying elliptic curves. In fact, it is impossible to obtain real points on the curves 
since we can't distinguish between a point and its opposite and we don't even know 
if we are on the curves or on their twists. If the elliptic curves are in the Weierstrass 
form, we only need the coordinates (x :: z) to test if the point is zero: just compute 
gcd (z,n). The morphisms between the Kummer surface and the underlying elliptic 
curves (identified with their quadratic twists) are rational over Q. However their 
building blocks are defined over an extension field of Q: they use square roots 
that disappear in the global morphisms. Since we are not allowed to use square 
roots we must rewrite the global morphisms. As a result their expressions become 
complicated however, they are computed only once for every run of HECM, thus 
they arc negligible before the cost of computing [k]P. Note that we only need 
the (x :: z) coordinates which are invariant under the isomorphism from the curve 
(in the Weierstrass form) to its twist. See the appendix for a description of the 
different elementary blocks. The first stage of the HECM algorithm is summarized 
in algorithm 4. 

We now turn to what is called stage 2. The initial point for the arithmetic is 
Q = [k}P, thus wc don't have the benefit of a "good" initial point. Moreover, 
stage 2 needs the x-coordinate of many points on the elliptic curve in Weierstrass 
form [13], therefore if we use hyperelliptic curves we need to apply morphisms a 
lot and, in this case, their cost would not be negligible. In addition to that, if 
we want to use Brent-Suyama's extension, we can't use pseudo-additions but real 
additions. For all these reasons, it seems that hyperelliptic curves should not be 
used for stage 2. Instead, we can apply the ECM stage 2 to the two underlying 
elliptic curves. 
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Algorithm 4 HECM (stage 1) 

Input: the number n to be factor. The smoothness bound B\. 
Output: a factor p of n. 
1: Compute k = 1cm (1,2,..., Bi). 

2: Choose a random decomposable curve C over Z/nZ and a point P on its Kum- 

mer surface. 
3: Compute Q = [k]P. 

4: Map Q to the two underlying elliptic curves Si. 

5: Hope that Q = mod p for one £j (test whether gcd (z, n) ^ 1). 

6: Else go to 2. 



3.2. Torsion of the curves. To improve the probability of success, we can force 
the group order of the elliptic curves modulo primes p to be divisible by small 
numbers. Mazur's theorem states that the torsion group E tor (Q) of any elliptic 
curve over Q is isomorphic to one of the following groups: 



In our case the underlying elliptic curves have at least four 2-torsion points which 
means that Z/2Z x Z/2Z is a subgroup of the torsion group. 

In theory it should be possible to find decomposable hyperelliptic curves with 
underlying elliptic curves having larger torsion group. However this leads to more 
complicated equations for the parametrization. There is also a theoretical obstruc- 
tion to improving the torsion: to have a torsion point of order more than 2 means to 
be able to find a starting point on the Kummer surface corresponding to a point on 
Jac(C) and not on its twist Jac{C). This theoretical obstruction could be removed 
if the curve and its twist had the same torsion. However this increases the number 
of equations needed for the parametrization. 

Using the parametrization presented in Section 2, we have table 1 where letters 
indicate whether and where 4-torsion points exist: C for the curve and T for its 
quadratic twist. With the same probability, we work on the curve or on its twist. 
Suppose that the Legendre symbols of s 2 — u 2 , s 2 — 1, u 2 — 1 are independent 
which experimentally is a valid approximation. Then if p = 1 [4] we have a point 
of 4-torsion with probability 1/2, and if p = 3 [4], the probability is 3/4. Computer 
experiments show that for p = 3 [4] the power of two dividing the cardinality is, on 
average, 3.48 (instead of 3.5) and for p = 1 [4] it is 3.15 (instead of 3). 

Experimentally, the order of the curves is as likely to be smooth as a random 
integer about 1/7.75 in value (compared with 1/23.7 for Suyama's curves [13]). 
The parametrization s = f^fy, u = 2, v = 1 + provides better torsion (the power 
of two is on average 3.66) but there are few such curves with low parameters. 

3.3. Small parameters. For a chosen Kummer surface and a chosen initial point, 
the arithmetic on Kummer surfaces uses many multiplications by fixed parameters. 
Remember from section 1.3 that the cost of a scalar multiplication is 4M + 125+16(i 
per bit of the multiplier; i.e., 4 multiplications, 12 squares and 16 multiplications by 
constants. Suppose that these constants are small in comparison to the number n to 




1 < to < 10 or m = 12 
1 < to < 4 
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Table 1 . 4-torsion points on the underlying elliptic curves 
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s 2 - u 2 


3 2 -l 


u 2 -l 


-1 


□ 


M 


□ 


□ 


□ 


c 


C T 




c 


C T 




□ 


c 






T 


C T 




□ 


□ 


T 


C T 




C 






□ 


C 


C T 




C 


C T 



be factored, then the cost of multiplications by them will be negligible with respect 
to the cost of the full-length multiplications. The cost per bit of the multiplier 
becomes 4M +125 and if we use the rule of thumb 15 = 0.8M then it is 13.6M. For 
one run, GMP-ECM (which uses Montgomery coordinates and the PRAC algorithm 
[7]) uses approximately 6M + 35. Thus two runs of ECM cost 12M + 65 or 16.8M 
with 5 = 0.8M. This is a speed-up of 20%. Of course this is only valid for large n. 
EECM (ECM with Edwards curves) uses signed sliding window to compute [k] P 
(see [1]). These chains use 1 elliptic curve doubling and e elliptic curve additions 
for each bit of k where e converges to when k increases. Doubling in Edwards 
coordinates uses only 3M + 45, while additions use 10M + 15 + Id. In theory, 
two runs of EECM would be faster than one run of HECM if e is less than 1/12 (if 
15 = 0.8M then e should be less than 1/20). However this requires that the width 
of the window be very large and thus precomputations and memory usage might be 
not negligible. The authors of [1] claim that with B x = 16384, EECM uses 195111 
multiplications in stage 1 (with a width-6 signed sliding window). For the same 
Bi, HECM uses 379334 multiplications but since it does two curves in parallel, it 
uses only 189667 multiplications for stage 1 for one elliptic curve. Therefore, more 
experimentations are needed but to our knwolegde there is no public implementation 
of EECM. 



By a small constant we mean a number which fits into a signed long. This con- 
straint limits the number of curves we can use: if we work with a 64 bits processor 
there are 185, 399 useful hypcrclliptic curves which is sufficient for finding factors 
of more than 65 digits. In the multiplication algorithm, the 16d could be reduced 
to 12<i since we work in the projective space but we would have rational numbers 
on Q that, modulo n, must each fit into a long which is more complicated. 
More explicitly, our parametrization gives (a : /3 : 7 : S) in terms of rational func- 
tions of small degree in (s,u) E Q 2 . The same is true for (1/A : l/B : l/C : 1/D). 
We work in projective coordinates so we can clear the denominators. The constants 
used in the multiplications become: 



'1111 

(3 7 5 



= (sV : s 5 v 2 : s (u 2 - s 2 ) (u 2 - l) : sV) 
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(l444) = ((-l)V-2^ 4 )( S + . 2 ) 2 : 

(u 2 -l) (s + u 2 ) 2 (s-u 2 ) 2 : 

-(« 2 -l) (s + u 2 f (s-u 2 ) 2 : 

(,s + l) 2 ( S 2 -2 U 2 + U 4 ) ( S -« 2 ) 2 ) . 

Note that s is rational so we take s = a/b for integers a, b. After generating a 
hyperelliptic curve, we test if these constants are small. We then need to find a 
point on the Kummer surface for which the inverses of its coordinates are small. 
Unfortunately we have found no point such that the polynomial expressions for 
the inverses of the coefficients have degree equal or lower than the degree of the 
constants above. In practice, we have several generic points so we check whether 
those points fit, if not we try with another curve. 

3.4. A numerical example. Let's try to factor n = 4,816,415,081 with HECM. 
We use 51 = 25 and B2 = 200. Let's use the following parameters: s — |, u = 2 
and v — 9. The Kummer surface is given by the parameters 

(a : (3 : 7 : 8) = (l : 2 : A : l) . 

One point on the Kummer surface is 

P = (-272 : 272 : 63 : -140). 

We compute [k] P on the Kummer surface with k = 1cm (2, Bl) = 26, 771, 144, 400 
and send the resulting point on the two underlying elliptic curves: 

(3455587574, 1) on 734346861 y 2 = (x - 1) (x ~ ~ > 

(3222355131, 1) on 2791313056 y 2 = (x - 1) (x - 25) (a; - 121) . 
We now use the ECM stage 2 on the two elliptic curves. The first one does not 
produce any factor but the second one produces the number 83003 which is a factor 
of n = 4816415081 = 58027 * 83003. 

We inspected the calculus and found that we didn't work on the hyperelliptic curve 
but on its twist: 

XV 2 =x{x-l)[x-^ (x-^j 

where \ is a non quadratic residue. The order of the initial point in the Jacobian 
of this curve modulo 83003 is 2 • 3 • 5 • 11 • 19 • 73 • 631. Its order on the first elliptic 
curve is 2 • 7 • 631 and on the second 2 • 3 • 5 • 19 • 73. 

4. Implementations and results 

There are many implementations of ECM: GMP-ECM is a free program based 
on the GNU MP library. It is described in detail in [13]. By using many functions 
from GMP-ECM, we have built a new program called GMP-HECM which uses hy- 
perelliptic curves. GMP-HECM will be distributed with GMP-ECM (currently it is 
distributed with the development version at http://gforge.inria.fr/projects/ecm/). 
Our software does stage 1 by generating a genus 2 hyperelliptic curve with small 
parameters, computes [k] P on the Kummer surface and send the resulting point 
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Table 2 . Optimal choice of the parameters for ECM and HECM 



Size 


GMP-ECM 


GMP-HECM 


of the 


Optimal 


Expected number 


Optimal 


Default 


Expected number 


factor 


Si 


of elliptic curves 


Bt 


B 2 


of elliptic curves 


1Q 2U 


11,000 


74 


14,000 


2.10 b 


75 


1Q 2b 


50,000 


214 


60,000 


16.10 b 


214 


1Q 3U 


250,000 


430 


260,000 


130.10 b 


491 


1Q 3b 


1.10 b 


904 


1.10 b 


900.10 b 


1,116 


1Q 4U 


3.10 b 


2,350 


3.10 b 


4.10 9 


2,871 


10 46 


11.10 b 


4,480 


11.10 b 


28. 10 9 


5,425 


10 50 


43. 10 b 


7,553 


43.10 b 


200. 10 9 


9,003 


10 55 


110.10 b 


17,769 


110. 10 b 


750. 10 9 


21,183 


10 bu 


260.10 b 


42,017 


260. 10 b 


2.10 ia 


49,534 


10 bb 


850.10 b 


69,408 


850. 10 b 


14.10 12 


81,387 



on the two underlying elliptic curves. Then, for stage 2, it uses GMP-ECM for the 
two underlying elliptic curves (This could be done in parallel). Indeed, GMP-ECM 
can use the results of stage 1 from an other software: it just needs the parameter A 
of the elliptic curve in short Weierstrass form y 2 = x 3 + Ax + B and the coordinates 
of the point on it. 

Our elliptic curves arc not, in general, Suyama curves so they don't have the 
same torsion groups as it is expected in GMP-ECM. This affects the probability 
of the cardinality of the curves being smooth and thus the choice of B\ and B 2 . 
GMP-ECM is able to choose the parameters B 2 for stage 2 from the value of B\ . 
This choice is not the same as the "standard continuation" (i.e. B 2 — 100 * B{) 
but produces much bigger B 2 . This reduces the number of expected curves for 
finding a factor. The choice of B\ is heuristic, we tried to minimize the means of 
the expected time. Note that for finding small factors, HECM is not interesting due 
to the cost of initialization and morphisms. Table 2 compared the optimal value 
of B\ for different sizes of the (usually unknown) factor. We see that for finding 
factors of at least 35 digits, we use the same B\ than ECM. Moreover the ratio of 
the expected number of curves needed for finding a factor decreases with the size 
of the factor. 

For not too small B\ , the theoretical and the real cost of the algorithm is linear 
in B\ . However its dependence in the size of the input number is not simple: modu- 
lar multiplications are quadratic for small n but become quasi-linear for large n. 
Moreover, there are a lot of other operations which complexity is not negligible for 
small n. Table 3 presents the time taken by the different operations during the 
computation of [k]P for different sizes of n. We see that squares take the same 
time as multiplications. The reason is that GMP-ECM has special assembly code 
for general modular multiplications but not for squaring and thus squarings use 
the same assembly code as normal multiplications. For multiplications by small 
constants, we made special assembly code. The table also shows that additions are 
not negligible for small n. 
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Table 3. Fraction of time taken by the arithmetic operations for 
different sizes of n 



size of n 


10 ioo 


10 lbu 


10 2UU 


1Q 2bU 


10 3U0 


10 350 


1Q 4U0 


10 bou 


10 iooo 


M 


0.143 


0.159 


0.194 


0.216 


0.209 


0.219 


0.216 


0.222 


0.234 


S 


0.428 


0.465 


0.557 


0.597 


0.622 


0.645 


0.65 


0.682 


0.700 


d 


0.163 


0.127 


0.092 


0.081 


0.070 


0.057 


0.037 


0.031 


0.051 


additions 


0.237 


0.179 


0.116 


0.109 


0.074 


0.057 


0.059 


0.041 


0.017 


S/M 


1.00 


0.98 


0.96 


0.93 


0.99 


0.98 


1.04 


1.02 


1.00 


d/M 


0.26 


0.20 


0.12 


0.093 


0.084 


0.066 


0.043 


0.035 


0.055 


Table 4. Comparison between stage 1 of GMP-ECM and GMP- 
HECM for different sizes of n with a constant B\ = 10 7 on a core 
2 at 2.4Ghz. 


size of n 


10 ioo 


10 ibu 


10 2.,u 


10 2b,) 


10 3UU 


1Q 3bU 


1Q 4UU 


10 buu 




HECM 
2*ECM 


1.15 


1.09 


1.03 


1.0 


0.97 


0.94 


0.95 


0.93 


0.89 



Table 4 compares the running time of two runs of GMP-ECM versus GMP- 
HECM (remember that one run of HECM is equivalent to two runs of ECM) for a 
fixed B\ and different sizes of n. This shows that for large n (at least 10 250 ), our 
software GMP-HECM is faster than GMP-ECM. Since HECM uses more squarings 
than ECM, having optimized assembly squaring code would be likely to provide a 
small improvement to the limit where GMP-HECM is more interesting than GMP- 
ECM. For very large n the speed up of GMP-HECM compared to GMP-ECM is 
around 11%. This result agrees with the theoretically predicted value. 

Conclusion 

Using a subfamily of hypcrclliptic curves of genus 2 (i.e decomposable curves), 
we have built an algorithm for factorization which is equivalent to two simultaneous 
runs of ECM. Kummer Surfaces allowed us to choose small parameters that give 
an efficient arithmetic in the Jacobian of the curve. In practice, for large numbers, 
our implementation is faster than GMP-ECM. Special assembly code for squaring 
would provide another improvement to GMP-HECM. 

Many thanks to Emmanuel Thome for his help on this work and for his comments 
on the draft versions. Thanks also to the developers of GMP-ECM and in particular 
to Alexander Kruppa for his great assembly code. 

Appendix: Morphisms 

A Kummer surface is defined by the first four theta constants (a : (3 : 7 : 5) . Its 
equation is 

4E' 2 af3~f6XYZT = ({X 2 + Y 2 + Z 2 + T 2 ) - 

-F {XT + YZ) - G (XZ + YT) - H (XY + ZT)f 
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where 



A = a + + i + S, 
C = a- + 1-6, 



E' 
F 
G 
H 



Y' = a + 0-i-6, 
T' = a- 0-1 + 6, 



ABCD 



{ad 


- 0i) {on - 


06) [a0 


{c? 


-0 2 -i 2 + 






{ad - 0i) 


5 


(a 2 


-0 2 +1 2 - 






{ai - 06) 


1 


{a 2 


+ /3 2 - 7 2 - 




(a0 - j6) 



In our case (decomposable curves with a = 5), the other theta constants are 
given by the following relations: 

0?(O) = a, 6l {O) = 0, 9 2 3 {0)=i, 6\{Q) = 5, 



9| (0) = e = ^Ja0-i5\ = Ja0 - i~6 

V a 



??o (0) 



ea a0 — j6 
Uf3 = e ' 



2 (0) = v 7 ^, 
CI2 



% (0) 



^(0) = 



?6 2 (0) 



75 — dip 

~oJW' 

ad — 0i 

~oJW' 



The global morphism in HECM goes from the Kummer surface to the product 
of the underlying elliptic curves modulo the isomorphisms which identify a curve 
with its twist. However we only need the (x :: z) coordinates which are invariant 
under this morphisms. It is the composition of elementary blocks which use square 
roots: 




C Jac{C) » Jac{C) / {±1} —^.JC 

f 

£\ x £ 2 



f. 

£\ x So- 



Let P = {X :Y:Z:T)& point on the Kummer surface IC^ a .p. T s) corresponding 
to the hyperelliptic curve C of equation y 2 — f (x) . We want the Mumford coordi- 
nates {u, v) of the two opposite divisors ip^ 1 (P) = {±D}. The divisors D and —D 
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share the same polynomial u and have opposite v polynomial. First, compute the 
following quantities where for readability we note 6% for Oi (0): 



2 7 (Z) 


-xelt 


q + Yeiel -zelel + T6l6l 






el (z) 


Afc, 6 t 'l0 - Y y 6 y 8 + Zy 5 y 8 ~~ J ^lO 






Oli (z) 




\rn2a2 _i_ yn2n2 rpn2n2 

- y v 6 v 10 + £v 5 v 10 - i o 5 v s 




ot-n 


*12 (Z) 








^8 ^ ^10 


0f 3 (z) 




?9 + yOjOfo + Z6pl - T9^6 2 Q 






2 u (z) 


-xen 


?| _ y^2^2 + ^2^2 + T6 2 d 2 






e\, (z) 


-xelt 


l 2 -Y0 2 6 2 + Ze 2 6 2 +T6 2 2 







In the general case, 8 2 e (z) is non zero. In this case u (£) = £ 2 + ui£ + uq is of 
degree 2 and ?;(£)= ± + v ) is of degree 1 with 



u = A 



2a2 (z) 

<M 6 (*)' 



y 8 u 14 1 



ui = (A - 1) 



gj_g?3 (£) 
<M 8 (*) 



- ?i - 1, 



#1#3#8 #14 ( z ) 



(e 2 e 2 8 2 e 2 16 ( Z )) 

+W 2 Q\Q\Q\ (XZ + YT) 



io£et e 2 ( Z ) e{ 2 ( z ) + eleje* e 2 (z) e{ x (z) 
e 2 e 2 



2 u 2n> n 2 , v , . . -t, 1 6 3 - w 2 w 1 , , ■) 



E> 



Z 2 + T 2 ) - 



-F {XT + YZ) - G (XZ + YT) - H (XY + ZT))) . 



Since u divides v 2 — / we obtain formulae to compute v\. If 8 2 6 (z) is zero and 
(A — 1) #5 9 2 3 (z) — A0f 9 2 4 (z) is non zero, then u (£) = £ + u a is of degree 1 and 
v (£) = ±v is constant with 



U 



\e 2 e 2 u {z) 



(\-i)e 2 e 2 3 (z)-xe 2 e 2 i (z) 



w(0 = ±«o = ±V/(-«o). 



The last case is when # 2 6 (z) = (A - 1) 6> 2 6> 2 3 (z) - A0f # 2 4 (z) = 0. In that case, the 
divisor D is the zero divisor. 



Let P be a point on the (2, 2)-decomposable hyperelliptic curve C given by the 
equation 

C : XV 2 = x(x — l) (x — A) (x — n) (x — v) , with A = /x- 



We want to map P to the elliptic curves. Let C be the curve given by 
C : «y 2 = (.t 2 - 1) (a: 2 - x\) (x 2 - x 2 3 ) 
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where 

/—, 7 A* + <7 1 - M - Q . . 

q = ±VM(/"-^), ^2 = — -, £3 = 1 + K = -X9/"(M- !)• 

The curves C and C are isomorphic by the change of coordinates 
C — » C 

fotf) >— » f^T^ W with W=7 ^ r 

The curve C maps to the elliptic curve £ : y 2 = (x — 1) (x — x 2 ) {x — xf) by the 
morphism (x,y) \— > (x 2 ,y). 

Changing g to — g changes the elliptic curve £ to the other underlying elliptic curve. 

Let / be the map from C to the product of the two elliptic curves £\ x £ 2 given by 
the product of the two maps defined above. The push forward /» of / is defined by 

, J Jac{C) — ► Jac{£\) x Jac(£ 2 ) 
U \ D = E[ =1 Pi r Poo ^ E[=i / (i'i) " r / (Poo) 

where we define / (Poo) = (Cfi , Cf 2 ) to be the zero of the two elliptic curves. Note 
that the divisors in the Jacobian of the elliptic curves are not reduced. For elliptic 
curves, the Jacobian of the curve is isomorphic to the set of points on the curve. 
Thus we can identify Jac(£i) with £j. We can rewrite /* as 

J Jac (C) — ► £\ x £ 2 
M £> = £•=! Pi -r^oo — E[=i/(P) 

A generic divisor D in the Jacobian of a genus 2 curve is the formal sum: D = 
Pi + P 2 — 2Poo. Thus in general we have to add the two points / (Pi) and / (P 2 ) 
on the elliptic curves. In practice, to do this operation, we need complete formulae 
for adding points on elliptic curves because we can't take a square root to get the 
y-coordinate which is needed for doubling in Wcierstrass coordinates. Thus we use 
Jacobi or Edwards equations for the curves. Another solution to avoid square roots 
is to work in an extension "field" of degree 4. 

For stage 2, GMP-ECM needs a point (x 2 ,y 2 ) on the curve y 2 — x 3 + Ax + B 
but we have a point (x\ :: z\) on nzy 2 = x 3 + a 2 x 2 z + a^xz 2 + a§z 3 . First translate 
the point by x ^ x — a 2 z/3 and divide the x-coordinate by z to get a point ?) 
on the curve ny 2 = f (x) = x 3 + a' 4 x + a' 6 . Then the points (x[, ±1) are on the 
curve 

Dy 2 = f{x) = x 3 + o! A x + a' e with D = f (x[) . 
By the change of variable 

tow) 

we get a point on the curve y 2 — x 3 + Ax + B with A = a' 4 / D 2 and B = ae/D 3 . 
Note that we can choose the sign of y since a point and its opposite have the same 
order (remember that we cleared the power of 2 in stage 1). 
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